/**************************************************************************//**
 * @file    biparabolic_interpolation.c
 * @author  İsmail SEZEN (email: sezenismail@gmail.com)
 * @author  Orlando Camargo Rodríguez (email: orodrig@ualg.pt)
 * @date    02-08-2009
 *
 * @note This is part of siklon project. (Siklon.c)
 *
 * <b>Description</b>
 *
 * Includes the routines to calculate the height at (x,y) point specified in a 3x3 matrix.\n
 *
 * Copyright (C) 2009 İsmail SEZEN
 *
 * <b>License</b>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *****************************************************************************/

/* Function Prototypes *///{
double get_product_for_val_base(double *x, double *var_i, int i_val);
double get_product_for_val(double *x, double *y, double *xj, int j_val,double *yi, int i_val);
double get_product_base(double *x,int i_val);
double get_product(double *xj,int j_val,double *yi,int i_val);
//}

/*! @brief Finds the function value of a specified x,y location by biparabolic interpolation
 *
 * @param[in] x is coordinate location on the x-axis
 * @param[in] y is coordinate location on y-axis
 * @param[in] *xj is a pointer to the array of x location values over the grid points
 * @param[in] *yi is a pointer to the array of y location values over the grid points
 * @param[in] zij is the 2D array represents the function values over the grid points represented by xj and yi arrays
 *  @return calculation result of the expression */
double Find_Zxy_Biparabolic_Interpolation(double x, double y, double *xj, double *yi, double *zij){
	register int i,j;//Loop variables
	double Zxy =0.0;//result value
	for(i=0;i<3;i++) {
		//printf("\n");
		for(j=0;j<3;j++) {
			double Zij = *( (double *) zij + ( (i*3) + j ) ); //Zij = Z( x[j], y[i] );
			//printf("[%g],",Zij);
			double aij = Zij / get_product(xj,j,yi,i); //aij is a coefficient for Zij
			double val = aij * get_product_for_val(&x,&y,xj,j,yi,i);//intermadiate value for Zxy addition process
			Zxy+= val;
		}
	}
	return Zxy;
}

/*! @brief  Subtracts each value of the 'xj' array (except 'j_val' element) from *x and multiplies them.\n
 *          Subtracts each value of the 'yi' array (except 'i_val' element) from *y and multiplies them.\n
 *          Multipiles all of the results obtained above.
 *
 *  For Example:\n
 *  for \f$ j_val = 1 \f$ and \f$ i_val = 2 \f$ \n
 *  @code result = (x - xj[0]) * (x - xj[2]) * (y - yi[0]) * (y - yi[1]); @endcode
 *
 *  @param[in]  *x          is pointer to the double precision value
 *  @param[in]  *y          is pointer to the double precision value
 *  @param[in]  *xj         is a pointer to the array
 *  @param[in]  *yi         is a pointer to the array
 *  @param[in]  j_val       is the reference point for the xj array
 *  @param[in]  i_val       is the reference point for the yi array
 *  @return calculation result of the expression */
double get_product_for_val(double *x, double *y, double *xj, int j_val,double *yi, int i_val){
	return get_product_for_val_base(x,xj,j_val) * get_product_for_val_base(y, yi,i_val);
}

/*! @brief  Subtracts each value of the 'p_arr' array (except 'landmark' element) from value and multiplies them.
 *
 *  For Example:\n
 *  for \f$ landmark = 1 \f$ ,\n
 *  @code result = (value - p_arr[0]) * (value - p_arr[2]); @endcode
 *
 *  @param[in] *value is pointer to the double precision value (x or y value)
 *  @param[in] *p_arr is a pointer to the xj or yi array
 *  @param[in] landmark is the reference point for the calculation
 *  @return calculation result of the expression */
double get_product_for_val_base(double *value, double *p_arr, int landmark){
	register int i;
	double result = 1.0;
	for(i = 0;i<3;i++)
		if( (i != landmark) ) result *= ( *value - *(p_arr+i) );
	return result;
}

/*! @brief A multiplication method
 *
 *  Subtracts each value of the 'xj' array (except 'j_val' element) from xj[j_val] element of the array and multiplies them.\n
 *  Subtracts each value of the 'yi' array (except 'i_val' element) from yi[i_val] element of the array and multiplies them.\n
 *  Multipiles all of the results obtained above.\n
 *  For Example:\n
 *  for \f$ j_val = 1 \f$ and \f$ i_val = 2 \f$ ,\n
 *  @code result = (xj[1] - xj[0]) * (xj[1] - xj[2]) * (yi[2] - yi[0]) * (yi[2] - yi[1]); @endcode
 *
 *  @param[in] *xj is a pointer to the array
 *  @param[in] *yi is a pointer to the array
 *  @param[in] j_val is the reference point for the xj array
 *  @param[in] i_val is the reference point for the yi array
 *  @return calculation result of the expression */
double get_product(double *xj,int j_val,double *yi,int i_val){
	return get_product_base(xj,j_val) * get_product_base(yi,i_val);
}

/*! @brief Subtracts each value of the 'p_arr' array (except 'landmark' element) from p_arr[landmark] element of the array and multiplies them.
 *
 *  For Example:\n
 *  for \f$ landmark = 1 \f$ ,\n
 *  @code result = (p_arr[1] - p_arr[0]) * (p_arr[1] - p_arr[2]); @endcode
 *
 *  @param[in] *p_arr   is a pointer to the array
 *  @param[in] landmark is the reference point for the p_arr array
 *  @return calculation result of the expression */
double get_product_base(double *p_arr,int landmark){
	register int i;
	double result = 1.0;
	for(i = 0;i<3;i++)
		if( (i != landmark) ) result *= ( *(p_arr+landmark) - *(p_arr+i) );
	return result;
}
